Stirling's approximations for exchangeable Gibbs weights* 



Annalisa Cerquetti 1 ^ 

Department of Methods and Models for Economics, Territory and Finance 
Sapienza University of Rome, Italy 

June 29, 2012 



Abstract 

We obtain some approximation results for the weights appearing in the exchange- 
able partition probability function identifying Gibbs partition models of parameter 
a £ (0, 1), as introduced in Gnedin and Pitman (2006). We rely on approximation 
results for central and non-central generalized Stirling numbers and on known re- 
sults for conditional and unconditional a diversity. We provide an application to an 
approximate Bayesian nonparametric estimation of discovery probability in species 
sampling problems under normalized inverse Gaussian priors. 

1 Introduction 

Exchangeable Gibbs partitions (Gnedin and Pitman, 2006) are the largest class of infinite 
exchangeable partitions of the positive integers N and are characterized by a consistent 
sequence of exchangeable partition probability functions (EPPF) in the form 



p(m,...,n k ) = V nik JJ(1 



for each n > 1, (m, . . . , n^) the sizes of the blocks, a £ (— oo, 1) and V n i~ coefficients 
satisfying the backward recursive relation 

V n ,k = (n- ka)V n+1)k + K+i,fe+i, 

for V\i = 1. The specific form of the V n & identifies the specific a-Gibbs model, each aris- 
ing (cfr. Th. 12, Gnedin and Pitman, 2006) as a mixture of extreme partitions, namely: 
Fisher's (1943) (a, £) partitions, £ = 1, 2, . . ., for a < 0, Ewens' (1972) Poisson-Dirichlet 
(9) partitions for a = 0, and Pitman's (2003) Poisson-Kingman PK(p a \t) conditional 
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partitions for a 6 (0, 1), where p a is the Levy density of the Stable subordinator. For 
a £ (0, 1), 9 > —a and mixing density 



70,a(*) 



r(e/a + i) M; ' 



where is the density of the stable distribution, the Poisson-Kingman (p a ,7a,e) 

model corresponds to the two-parameter Poisson-Dirichlet partition model (Pitman, 
1995, Pitman and Yor, 1997) with EPPF in the form 



p a> e(ni, ...,n k ) 



(9 + a)fc-ita 

(0 + l)n-l 



which reduces to the Dirichlet (9) partition model for a = 0. Recently exchangeable 
Gibbs partitions and the corresponding discrete random probability measures, have 
found application as prior models in the Bayesian nonparametric treatment of species 
sampling problems, (cfr. Lijoi et al. 2007b, Lijoi et al. 2008). Here interest typically lies 
in estimating the diversity of a population of species with unknown relative abundances, 
by estimating both predictive species richness (cfr. Favaro et al. 2009) and posterior 
species evenness (cfr. Cerquetti, 2012). By assuming the sequence of labels of different 
observed species to be a realization of an exchangeable sequence (Xj)j>i with de Finetti 
measure belonging to the Gibbs class, posterior predictive results are obtained with re- 
spect to a future m-sample conditioning on the multiplicities (ni, . . . ,n^) of the first k 
species observed in a basic n-sample. 

With respect to the whole PK(a, 7) class, the (a, 9) partition model ([I]) stands out 
for its mathematical tractability, and a huge amount of results are available for this 
model both in random partitions theory like in Bayesian nonparametric implementa- 
tions. But in general, without loss of generality, for mixing distributions in the form 
j(t) = h(t)f a (t), for f a (t) the Stable density, the Gibbs weights are obtained by mix- 
ing the conditional EPPF p a (ni, . . . ,nk\t), (cfr Pitman, 2003, eq. (66)) by the specific 
mixing distribution, which results in the following integral 



V 



a,h 
n.k 



a 



-ka 



T(n — ka) 



P 



n—l—ka 



/«((! - p)t)dp 



h(t)dt. 



(2) 



It follows that calculating explicitly the PK(p a ,h x f a ) coefficients outside the {a, 9) 
class can be hard, mostly due to the lack of explicit expression for the a-stable density. 
To give an example, in the Bayesian nonparametric implementation of those models, 
the most studied alternative to the (a, 9) model is the normalized generalized Gamma 
class (Pitman, 2003; see Lijoi et al. 2005, 2007a, Cerquetti, 2007) whose mixing density 
arises by the exponentialy tilting of the stable density, hence 



7 (t) = h(t) x f a (t) = exp{^(A) - Xt}f a (t), 
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for ip a (X) = (2A) a , A > the Laplace exponent of f a . By the reparametrization A = 

/3 1 /"/2, 

laAs~ 1/a ) = exp j/3 - \ 0) Us- l ' a )a- l s- l / a - 1 (3) 
and the corresponding iy n ^Y (Cerquetti, 2007) are given by 



r (") 7o + 2X) n ~ ka 

and can be rewritten (Lijoi et al. 2007a) by the change of variable x = (/3 1//q + 2A) Q , 

as a linear combination of incomplete Gamma functions 

v ; i=0 v 7 

The computational burden implicit in working with those priors is then clear. Addition- 
ally, when dealing with posterior predictive analysis in this setting, Bayesian nonpara- 
metric distributions for quantities of interest in species sampling problems are typically 
characterized by a ratio, that we term posterior Gibbs weights, of the kind 



V>n+m,k+k* 



V, 



n.k 



(5) 



for k* the number of different species observed in the additional m-sample, whose com- 
plexity increases with the complexity of the prior weights. It looks therefore interesting 
to investigate the possibility to obtain some sort of approximation result for both prior 
and posterior Gibbs coefficients. Here we propose a Stirling's approximation by relying 
on approximation results for both generalized central and non central Stirling numbers. 
The paper is organized as follows. In Section 2. we provide some preliminaries and 
obtain an approximation result for the prior Gibbs weights of the PK(p a , h x f a ) class. 
In Section 3 we obtain the approximation result for the posterior Gibbs weights (0) by 
relying on a recent result for conditional a diversity (Cerquetti, 2001) and in Section 4. 
we apply our findings to an approximate posterior estimation of discovery probability 
under normalized inverse gaussian partition model (Pitman, 2003, Lijoi et al. 2005) 
which corresponds to the generalized Gamma case for a = 1/2. 



2 Stirling's approximations for PK(p a , h x f a ) weights 

From the first order Stirling's approximation for ratio of Gamma functions r("+6) ~ n a ~ b 
is an easy task to verify that, for n — > oo and k ~ sn a , the following approximation 
holds for the EPPF ([I]) of the PD(a, 9) model 

. (af-^ik) IY0 + 1) fk\ e,a ^j r 

r(„) TWaA) na-**-'- W 



3 



Our first aim is to generalize ([6]) to the PK(p a , , y) class. First recall that general- 
ized Stirling numbers are combinatorial coefficients appearing in the expansion of ris- 
ing factorials {x) n = x(x + 1) ■ ■ ■ (x + n — 1) in terms of generalized rising factorials 
(x)k\a = x ( x + ot){x + 2a) ■ ■ ■ (x + (k — l)a), hence 



k=l 

An explicit expression for those numbers is given by Toscano's (1939) formula 

k 



S 



-1,-Q 

n,k 



a k k\^ ^ 



, (k 



-ja) 



(7) 



additionally (cfr. e.g. Pitman, 2006) they correspond to partial Bell polynomials of the 
kind 



B n ,k(l ~ «.) 



in 



{A u ...,Ak}€V [n]k j=l ( ni ,...,n k )j=l 



j2 -q( 1 -«)%~i 



for a £ (0, 1), where the first sum is over all partitions (A±, . . . , Ak) of [re] in k blocks, 
and the second sum in over all compositions of re into k parts. The following result relies 
on an approximation for generalized Stirling numbers obtained in Pitman (1999) and 
on a local approximation result for the distribution of a proper normalization of the 
numbers of blocks K n of an a £ (0, 1) Gibbs partition. 

Proposition 1. Let be the coefficients in the EPPF of an exchangeable Gibbs 

partition arising by a general PK(p a , h x f a ) Poisson-Kingman model, then the follow- 
ing Stirling's approximation holds for re — > oo and k ~ sn a 



V 



n,k 



a k - L T(k) 
T(n) 



k 



ii L 



-l/a 



Proof. In Pitman (1999, cfr. eq. (96)) an asymptotic formula for the generalized Stirling 
numbers S~^~ a for n — > oo and < s < oo , with k ~ sn a is derived by known local 
limit approximations for the number of blocks in a partition generated by a PD(a,a) 
model by a stable density, namely 

S n,k ~ r(fc) Sa{s)n , (8) 

where, g a (s) = a~ 1 f a (s~ 1 ^ a )s~ 1 ~ 1 ^ a . Now, for K n the number of blocks in a PK(p a , hx 
f a ) partition model (Pitman, 2003) almost surely, for n — > oo 

K n /n a -> S a 
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for S a ~ h(s 1 / a )g a (s), which implies the following local limit approximation for k ~ 
sn a holds for the distribution of K n 

F(K n = k) « V:£s^- a « 

Thus the result follows by substitution. □ 

Example 2. The result in Proposition 1. agrees with the Stirling's approxima- 
tion for the weights of the (a, 9) model obtained in Q. It is enough to notice that 
PD(a,8) = PK(p a ,'y ai g) for ^g jCt = h a> g x f a corresponding to the 6 polynomial tilting 
of the stable density for 

hit) = Tid + 1) t~° 

As for the Generalized Gamma model recalled in (|3j) and @, the first order Stirling's 
approximation for the Gibbs weights will correspond to 



T(n) 

3 Stirling's approximation for posterior PK(p a , hxf a ) weights 

As recalled in the Introduction, posterior predictive distributions in Bayesian nonpara- 
metric estimation in species sampling problems under Gibbs priors are typically charac- 
terized by a ratio of Gibbs weights. As an example consider the posterior joint distribu- 
tion of the random vector K m , L m , Si, ... , Sx m for K m the number of new species gen- 
erated by the additional m sample, L m the number of new observations in new species, 
and Si, ... , SK m the vector of the sizes of the new species in exchangeable random order, 
namely 

¥(K m = k,L m = s,S Km = (si,...,s Km )\n) = 

k* 



S\ Vn +m,k+k* I Tn> 

si\---s k *\k\ V n ,k \s y 

1=1 



J (n - ka) m - s - a) Si -i 



or the Bayesian nonparametric estimator for the number K m of new species in the 
additional sample (Lijoi et al. 2007b) 

m T/ 

E(K m \K n = k)=J2 k* Vn+ ™> k+k * s-\: a ^ n ~ ka \ 

k*=o Vn ' k 

Here 5' m 1 ^ a ' ^ n kc ^ are non-central generalized Stirling numbers, defined by the con- 
volution relation (see Hsu and Shiue, 1998, Eq. (16)) 

m / \ 

s- m ]k-* a >- {n - ka) = E (?)(»- k *)m-ss;^ a , (9) 



s=k* 
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whose corresponding Toscano's formula ([7]) may be derived by ([9]) as follows 



. s J a* k*l ^— ' V 7 



^ E(-D J ) E (T) (n-fea) m -,(-ja). = ~ Q (*-Cj+*)«), 

(10) 



3=1 s=fc* 3=1 



In order to derive an approximation result for 



we replicate the approach adopted in the previous Section by first looking for an ap- 
proximation for non central generalized Stirling numbers by exploiting their definition 
in terms of convolution of central generalized Stirling numbers © • 



Proposition 3. The following approximation holds for non central generalized Stir- 
ling numbers SLA* a ' for a € (0, 1) and m large 

q-l.-a-Cn-to) a 1 r(m) n _ ka _ a f , 1 s n _fc a _i -a-1 f-.-ay,. f 1 -i \ 

/or 

and z = k*/m a . 

Proof. By the definition of non-central Stirling numbers 

m / \ 

S:}C^- {n - ka) = ( m )(n-ka) m -.^- a 



and by 



This may be rewritten as 



9 -l -a ,-(n-to) _ i_fc* mr(m) r(g) T(n - fcq + m - s) / k* m c 

m > k * " ^ T(s)r(m - s + 1) T(k*) T{n-ka) S 9a \m a s a 

For z = k* /m a , by first order Stirling approximations for ratio of Gamma functions 

a-l,-a-(n-ka) _ a 1 (mjm / m - s \ (~ n -<*\(JL 

m > k * ~ r(n- fca)r(Jfe*) m / 



s=fc* 
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and by the change of variable p = s/m 



a 1 k * m a+n ka T(m 
T(n- ka)T(k*) 



A- 

Jo 



v) n ~ ka - x g a (,z V - a ) v - a - x d v . 



Now we are in a position to obtain the Stirling's approximation for the ratio V n +m,k+k* /Vn,k 
resorting to a general result for conditional a diversity for exchangeable Gibbs partitions 
driven by the Stable subordinator, recently obtained in Cerquetti (2011). 

Proposition 4. The following Stirling's approximation holds for the general posterior 
Gibbs weights for a PK(p a , h x f a ) partition model for large m and k* ps sm a : 

Vn+m,k+k* _ fr(s- 1/a )s fc a fc *r(fc*)r(n)m-( ra - fcQ ) _ a k+k *- 1 h(s- 1 / a )s k T(k*)m-( n - ka ') 
Vn,k ~ E nAa (h(S~V«))T(m)T(k) " V>) ' 

(12) 

Proof. Let II be a PK(p a ,~/) partition of N driven by the stable subordinator for some 
< a < 1 and some mixing probability distribution that without loss of generality we 
assume in the form jit) = h(t)f a (t). Fix n > 1 and a partition (m, . . . , n^) of n with 
k positive box-sizes, then by a result in Cerquetti (2011) for K m the number of new 
blocks induced by an additional m-sample 



for S™ ,k , having density 



for 



^\{K n = k)^sH 



Jn ' k[ > E« fc [/i(5-V«)] ' { ' 



the density of the product Y a ^ x [W] a where Y a ,k has density 

. , T(ka + 1) k . . 

g a ,ka\y) = Y(k + i) y 9a ^ y ' 

for g a (y) = ce~ 1 y~ 1 ~ 1 ^ a faiu^ 1 ^ 01 ), independently of W ~ /3(ka,n — ka). Additionally 

K,kM b >i ~ Vn > k < h — frfy — • 

Hence the following local approximation for the posterior distribution of K m (Lijoi et 
al. 2007) holds 

F ( K -h*\K - b\ - Vn + m ' k + k * q -l,-a,-(n-ka) _ h ( s lA *)ffn,fc( s ) a 



□ 
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Substituting (fTT|) the result easily follows. □ 

Example 5. [Poisson-Dirichlet (a, 9) model] Direct first order Stirling's approxima- 
tion for the posterior PD(a, 6) coefficients provides 

Vn +m ,k+k* _ (e + ka) k ^ a _ a k *({k*) e / a+k )T{k*)T(9 + n) 
V n , k {0 + n) m ~ T{m)T(9/a + k)m d + n ' 

while an application of (fl~2|) for h(s~ l / a ) = r(e/a+i) y^ds 



n+m,k+k* 



V, a 



(^)" Va J ^) S r(F)m>-Hr(fl + l) T(9/a + l)T(9 + 



,fc. 

II/,- I i n I i H 4- 1 i 

n) 



V n>k T{m)T(9/a + l) a k ~ l T(9/a + k)T{9 + 1) 

and the two results agree. 

Example 6. [Generalized Gamma (/?, a) model] The exact form of the ratio in this 
case is given by 

y n,k (n) m Er=o (V) (-iyw /a nk - p) 

while an application of (|12p with 

h[{k* /m*)- 1 /*] = exp (j^j 

provides 

a k* expj-f (|r) 1/Q | (^) fc r(fe>-("- to )r(„) 



V r 



n+m,k+k* 



V n , k r(m) E?=o (V) i-lWV-nk - ±; /J) 

4 Approximate estimation of discovery probability under 
inverse Gaussian partition model 

In Lijoi et al. (2007b) a Bayesian nonparametric estimate, under squared loss function 
of the probability of observing a new species at the (n + m + l)th draw condition- 
ally on a basic sample with observed multiplicities (ni, . . . ,n k ), without observing the 
intermediate m observations has been obtained as 



rSn.fc _ \ Vn+m+l,k+k*+l c -l,-a,-{n-ka) 
k*=0 Vn ' k 



The authors even derive explicit formulas under Dirichlet, two-parameter Poisson-Dirichlet 
and normalized Inverse Gaussian priors, which correspond to the generalized Gamma 
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prior model for a = 1/2. Notice that the larger is the size of the additional sample the 
more computationally hard would be to calculate explicitly (|14|) . 

By an application of (|12p . an approximate evaluation of (j!4[) for large m would 
corresponds to 

M „ V a k * +k Ks- l l a )s k nk* + l)(m+ i,- a ,- (w - fca) _ 

m ~ fe ^ K, fc r(m + 1) m > fc * 

rv^tm _l_ ~\\~ ( n— m 

= rrri-^n E fc((W 1/o )(WVr(f + i)^ ta) , 

which specializes under inverse Gaussian model as follows 

fyn,k _ (m+l)-("- fe / 2 )r(n) 

-IS r 



EJ 771+1 / /3\ 2 \ k* -l,-l/2,-(n-2k) 

and by the Toscano's formula for generalized Stirling numbers recalled in (|10|) simplifies 
to 

£n, fe = ( m+ l)-("-fc/2)r(n) 

m r(m + 1) Er=o ("7 1 ) (-i) l /? 2i r(fe - 2z ; /?) 

To have an idea of the reduction of the computational burden obtained provided by the 
approximation here we report the exact formula obtained in Lijoi et al. (2007b) 

m = (-P 2 r +i a Er= m rr) m 2 )^ + k* + 1 + 2 * - 2 (m + n) ; # w 



,s J \ s — I / !'(/,' 



/ m\ /2s - fc - 1\ 2 fc " 2s r(s) 

EL «_i Hy^"-^ 



s=fc 
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